%%%%%%%%%%%%%%%
% Housekeeping
%%%%%%%%%%%%%%%

clc;
clear                                   all;
close                                   all;
beep                                    off;
rng.Seed                                = 110887;
set(groot,'defaultAxesFontName', 'Garamond')
set(groot,'defaultTextFontName', 'Garamond')
set(groot,'defaultAxesTickLabelInterpreter','latex'); 
set(groot,'defaultLegendInterpreter','latex');
cd(pwd); 
addpath(strcat(pwd,'/auxiliary_functions/'));
project_folder                          = fileparts(fileparts(pwd));
root_folder                             = fileparts(project_folder );

%%%%%%%%%%%%%%%
% Model setup
%%%%%%%%%%%%%%%

% Structural parameters
pars.r                                  = -1/12*log(0.70);                                                              % Monthly model with annual discount rate of 0.70 
pars.k                                  = 0.163;                                                                        % Upper bound on adjustment intensity
pars.M_c                                = 1;                                                                            % Cash stock average
pars.M_e                                = 0.97;                                                                         % Electronic money    
pars.C                                  = 0.060;                                                                        % Strength of complementarities           
pars.mean_reversion                     = 120;                                                                          % Shock persistence, expressed as the number of days to get back to within 10% of initial value
pars.theta                              = -log(1-0.90)/( pars.mean_reversion/30 );                                      % Shock persistence, expessed as OU parameter
pars.sigma                              = 0.060;                                                                        % Standard deviation of money shock    
pars.T                                  = 12*10;                                                                        % Date after which mean-reversion goes to zero
pars.tIRF                               = 12*2;                                                                         % IRF horizon

% Discretization parameters
Nu.NX                                   = 1e2+1;                                                                        % # of gridpoints for X
Nu.N                                    = 1e2;                                                                          % (1/2)*( # of gridpoints for M - 1 )
Nu.Nmin                                 = ceil(pars.theta/pars.sigma^2* ...
                                               ((pars.r+pars.k+pars.theta)/(pars.r+pars.k))^2* ...
                                              max((pars.M_c-pars.M_e)^2,(pars.M_e + pars.C - pars.M_c)^2)); 
Nu.N                                    = max(Nu.N,Nu.Nmin); 
Nu.Deltat                               = 1/(Nu.N*pars.theta);                                                          % Time period, expressed as fraction of a month
Nu.theta_target                         = pars.theta;
Nu.solutionmode                         = 1;

%%%%%%%%%%%%%%%%%%%%%
% Solution for C = 0
%%%%%%%%%%%%%%%%%%%%%

pars.C                                  = 0;
solution_noC                            = solve(pars,Nu,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Boundaries to illustrate
% the solution with fixed
% costs
%%%%%%%%%%%%%%%%%%%%%%%%%%%

M_bar                                   = pars.M_e - pars.theta/(pars.r+pars.k)*(pars.M_c - pars.M_e);
M_s                                     = M_bar + 0.02;
M_S                                     = M_bar + 0.18;                               

%%%%%%%%%%
% Plots
%%%%%%%%%%
fignum                                  = 1;
tvec                                    = (-solution_noC.preT.Nu.Deltat):solution_noC.preT.Nu.Deltat:solution_noC.preT.pars.T;
tvecIRF                                 = (-solution_noC.preT.Nu.Deltat):solution_noC.preT.Nu.Deltat:solution_noC.preT.pars.T;
color1                                  = [0, 0.4470, 0.7410];
color2                                  = [0.6350, 0.0780, 0.1840];
gray                                    = [0,0,0]+0.3;
gray2                                   = [0,0,0]+0.5;
gray3                                   = [0,0,0]+0.7;
color3                                  = [0.9290, 0.6940, 0.1250];
color4                                  = [0.4660, 0.6740, 0.1880];

% Phase diagram - Large shock
ff2 = figure(fignum);
fignum = fignum+1;

    X                                   = solution_noC.preT.Nu.X';
    Phi                                 = solution_noC.preT.Phi_0;
    Phiref                              = Phi(1);
    ylim1                               = 0.8;
    ylim2                               = 1.2;
    Deltat                              = solution_noC.preT.Nu.Deltat;
    k                                   = solution_noC.preT.pars.k;
    theta                               = solution_noC.preT.pars.theta;
    M_c                                 = solution_noC.preT.pars.M_c;
    NX_q                                = 1e1;
    NM_q                                = 1e1;
    Xgrid                               = repmat(linspace(0,1,NX_q),[NM_q 1]);
    Mgrid                               = repmat(linspace(ylim1,ylim2,NM_q)',[1 NX_q]);
    Phi                                 = solution_noC.preT.Phi_0;
    Phigrid                             = interp1(X,Phi,Xgrid(:));
    Phigrid                             = reshape(Phigrid,size(Xgrid));
    dXgrid                              = k*Deltat*( (1-Xgrid).*( Mgrid <= Phigrid ) - Xgrid.*( Mgrid > Phigrid ) );
    dMgrid                              = theta*Deltat*( M_c - Mgrid );
    index_above                         = Mgrid > Phigrid;
    index_below                         = Mgrid <= Phigrid;

    X                                   = solution_noC.preT.Nu.X';
    Phi                                 = solution_noC.preT.Phi_0;
    Phiref                              = Phi(1);
    tPlot                               = solution_noC.IRF.tvec;
    ylim1                               = 0.8;
    ylim2                               = 1.2;
    M_0                                 = pars.M_c;
    M_t                                 = pars.M_c*(1-exp(-pars.theta*tPlot)) + exp(-pars.theta*tPlot)*(1-0.175)*pars.M_c;
    X_0                                 = 0;
    X_t                                 = exp( - (pars.k)*tPlot )*X_0 + (1-exp( - pars.k*tPlot ));
    indexswitch                         = find(M_t >= M_s,1,'first');
    X_t(indexswitch:end)                = X_t(indexswitch);
    
    plot(X,M_s*ones(size(X)),'-','color',gray2,'linewidth',0.5); 
    hold on;
    plot(X,M_S*ones(size(X)),'-','color',gray2,'linewidth',0.5); 
    plot(X,pars.M_c*ones(size(Phi)),'--','color',gray3,'linewidth',0.5);
    hold on;
    ylim([ylim1,ylim2]);
    xlabel('$X_t$','interpreter','latex');
    ylabel({'$M_t \quad$ ',''},'interpreter','latex','Rotation',0);
    th1 = text(1+0.04,pars.M_c,{'$\quad M^c$'},'interpreter','latex','HorizontalAlignment','center','color',gray3,...
                                'VerticalAlignment','middle','fontsize',12); 
    th2 = text(1+0.04,M_s,{'$\quad M_s$'},'color',gray2,'interpreter','latex','HorizontalAlignment','center',...
                                'VerticalAlignment','middle','fontsize',12); 
    th2 = text(1+0.04,M_S,{'$\quad M_S$'},'color',gray2,'interpreter','latex','HorizontalAlignment','center',...
                                'VerticalAlignment','middle','fontsize',12); 
    set(th2,'Rotation',0)
    ax = gca;
    ax.YLim = [ylim1 ylim2];
    ax.YTick = [];
    ax.XLim = [0.00 1.00];
    rna = fill([X' flip(X)'],[(M_S*ones(size(X)))' (ax.YLim(2)*ones(size(X)))'],color3,'LineStyle','none');
    rna.FaceAlpha = 0.15;
    rna = fill([X' flip(X)'],[(M_s*ones(size(X)))' (ax.YLim(1)*ones(size(X)))'],color4,'LineStyle','none');
    rna.FaceAlpha = 0.15;
    rna = fill([X' flip(X)'],[(M_s*ones(size(X)))' (M_S*ones(size(X)))'],gray2,'LineStyle','none');
    rna.FaceAlpha = 0.1;
    text(0.5,Phiref-0.025,{'$dX_t > 0$'},'color',color4,'interpreter','latex',...
                    'HorizontalAlignment','center','Fontsize',10);  
    text(0.5,Phiref+0.225,{'$dX_t < 0$'},'color',color3,'interpreter','latex',...
                    'HorizontalAlignment','center','Fontsize',10);
    text(0.5,Phiref+0.075,{'$dX_t = 0$'},'color',gray3,'interpreter','latex',...
                    'HorizontalAlignment','center','Fontsize',10); 

    plot(X_0,M_0,'o','MarkerSize',10,'MarkerFaceColor','white','MarkerEdgeColor',gray); hold on;
    plot(X_t(1),M_t(1),'o','MarkerSize',10,'MarkerFaceColor',gray,'MarkerEdgeColor',gray); hold on;
    set(findall(gcf,'-property','FontSize'),'FontSize',25)    

    set(ff2,'Units','inches');
    screenposition                              = get(ff2,'Position');
    set(ff2,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff2.PaperPositionMode                       ='auto';
    set(ff2,'position',[0,0,8,7])
        X_t                                         = X_t(1:(indexswitch-1));
        M_t_temp                                    = M_t(1:(indexswitch-1));
        p1 = plot(X_t,M_t_temp,'-','color',gray,'LineWidth',2,'MarkerSize',0.01,'MarkerFaceColor',gray2,'MarkerEdgeColor',color1); hold on;
        X_t                                         = exp( - (pars.k)*tPlot )*X_0 + (1-exp( - pars.k*tPlot ));
        X_t(indexswitch:end)                        = X_t(indexswitch);
        p2 = plot(X_t,M_t,'-','color',gray,'LineWidth',2,'MarkerSize',0.01,'MarkerFaceColor',gray2,'MarkerEdgeColor',color1); hold on;
        line2arrow(p1);
        line2arrow(p2);
    exportgraphics(ff2,[root_folder '/output/figures/' 'Figure_H_15.pdf'],'BackgroundColor','none');